% Table 5 in "Commodity Prices, Convenience Yields, and Inflation,"
% March 2011

clear all;

load irates_monthly_oecd.dat;   % other G-7 interest rates
load erates_monthly_oecd.dat;   % other G-7 exchange rates against USD
load cpi_monthly_oecd.dat;      % other G-7 CPI indeces (all items)

for country=1:6;
i=irates_monthly_oecd(:,country);
P=cpi_monthly_oecd(:,country);
er=erates_monthly_oecd(:,country);
nn=rows(P);

h=1;            % forecast horizon
rhat=2;         % number of principal components
nwl=-1;         % number of lags in Newey-West HAC estimator
bsize=4;        % bootstrap block size
B=4999;         % number of bootstrap replications

if country==1; load cyj_can.dat; cyj=cyj_can; load qj_can.dat; qj=qj_can;
elseif country==2; load cyj_jap.dat; cyj=cyj_jap; load qj_jap.dat; qj=qj_jap;
elseif country==3; load cyj_fra.dat; cyj=cyj_fra; load qj_fra.dat; qj=qj_fra;
elseif country==4; load cyj_ger.dat; cyj=cyj_ger; load qj_ger.dat; qj=qj_ger;
elseif country==5; load cyj_ita.dat; cyj=cyj_ita; load qj_ita.dat; qj=qj_ita;
elseif country==6; load cyj_uk.dat; cyj=cyj_uk; load qj_uk.dat; qj=qj_uk;
end

[ehat,Fhat,lamhat,ve]=pc(standard(cyj),rhat);       % PC estimation
[ehat2,Fhat2,lamhat2,ve2]=pc(standard(qj),rhat);    % PC estimation

infl_h=100*(log(P(3+h:nn,1))-log(P(3:nn-h,1)));
infl_1=100*(log(P(3:nn-h,1))-log(P(2:nn-h-1,1)));
infl_2=100*(log(P(2:nn-h-1,1))-log(P(1:nn-h-2,1)));
p1_cy=Fhat(3:nn-h,1);
p2_cy=Fhat(3:nn-h,2);
p1_rcpd=Fhat2(3:nn-h,1);
p2_rcpd=Fhat2(3:nn-h,2);
roil1=qj(3:nn-h,22);
roil2=qj(2:nn-h-1,22);
dx=100*(log(er(3:nn-h,1))-log(er(2:nn-h-1,1)));
ir=i(3:nn-h,1);

% Inflation equation
x1=[ones(rows(p1_cy),1) p1_cy p2_cy p1_rcpd p2_rcpd roil1 roil2 infl_1 infl_2];
res1=nwest(infl_h,x1,0);
b1=res1.beta;
r1=res1.resid;
R2_1=res1.rbar;
G1=nw((r1*ones(1,cols(x1))).*x1,nwl);
V1=inv(x1'*x1)*G1*inv(x1'*x1);
se1 = sqrt(diag(V1));
ts1 = b1./se1;

% Bootstrap procedure
bmatrix=[cyj(3:nn-h,:) qj(3:nn-h,:) roil1 roil2 infl_h infl_1 infl_2 p1_cy p2_cy p1_rcpd p2_rcpd];
g=floor(rows(bmatrix)./bsize);
bmatrix=bmatrix(1:g*bsize,:);
betab=[]; tstatb=[];
for j=1:B;
    indexx=floor(1+g*unifrnd(0,1,1,g))';
    ind=blockind(indexx,bsize);
    ind=ind(1:rows(bmatrix));
    matrixb=bmatrix(ind,:);
    cyjb=matrixb(:,1:23);
    qjb=matrixb(:,24:46);
    roilb1=matrixb(:,47);
    roilb2=matrixb(:,48);
    inflb_h=matrixb(:,49);
    inflb_1=matrixb(:,50);
    inflb_2=matrixb(:,51);
    p1_cyt=matrixb(:,52);
    p2_cyt=matrixb(:,53);
    p1_rcpdt=matrixb(:,54);
    p2_rcpdt=matrixb(:,55);
    
    [ehatb,Fhatb,lamhatb,veb]=pc(standard(cyjb),rhat);
    [ehatb2,Fhatb2,lamhatb2,veb2]=pc(standard(qjb),rhat);
    c1=corr([Fhatb(:,1) p1_cyt]);
    c2=corr([Fhatb(:,2) p2_cyt]);
    c3=corr([Fhatb2(:,1) p1_rcpdt]);
    c4=corr([Fhatb2(:,2) p2_rcpdt]);
    p1_cyb=sign(c1(1,2))*Fhatb(:,1);
    p2_cyb=sign(c2(1,2))*Fhatb(:,2);
    p1_rcpdb=sign(c3(1,2))*Fhatb2(:,1);
    p2_rcpdb=sign(c4(1,2))*Fhatb2(:,2);
    
    xb=[ones(rows(p1_cyb),1) p1_cyb p2_cyb p1_rcpdb p2_rcpdb roilb1 roilb2 inflb_1 inflb_2];
    [bb,bintb,rb,rintb,statsb] = regress(inflb_h,xb);
    Gb=nw((rb*ones(1,cols(xb))).*xb,nwl);
    Vb=inv(xb'*xb)*Gb*inv(xb'*xb);
    seb = sqrt(diag(Vb));
    betab(j,:)=(bb-b1)';
    tstatb(j,:)=((bb-b1)./seb)';
end;
fprintf(' h = %6.4f\n',h);
fprintf(' country = %6.4f\n',country);
fprintf('    coeff.     s.e.    boot s.e.  90%% bootstrap CI\n')
% equal-tailed percentile-t bootstrap
disp([b1 se1 std(betab)' b1-se1.*quantile(tstatb,0.95)' b1-se1.*quantile(tstatb,0.05)'])
% Hall's percentile bootstrap
%disp([b1 se1 std(betab)' b1-quantile(betab,0.95)' b1-quantile(betab,0.05)'])
fprintf(' Rbar^2 = %6.4f\n',R2_1);

end;